#!/usr/bin/python
# -*- coding: utf-8 -*-

import os
import struct
import array

import scipy
import numpy

import lib
import nmr_prototype
import peaks
import integrals

# Variablen, die man aus den Dateien suchen soll
# Im Format name : typ
searchAcqu = {
    'SFO1'    : numpy.float64,
    'SFO2'    : numpy.float64,
    'SFO3'    : numpy.float64,
    'SFO4'    : numpy.float64,
    'SFO5'    : numpy.float64,
    'SFO6'    : numpy.float64,
    'SFO7'    : numpy.float64,
    'SFO8'    : numpy.float64,
    'SW'      : numpy.float64, # Spektrumbreite (ppm)
    'SF'      : numpy.float64, # Im Prinzip sollte hier die Frequenz von 0ppm stehen.
    'SOLVENT' : str, # Das verwendete Lösemittel.
    'GB'      : numpy.float64, # GM-Faktor => Window-Funktion
    'LB'      : numpy.float64, # Linienverbreiterung (Hz) => Window-Funktion
    'TD'      : numpy.int32, # Acquirierte Punkte (r+i), also Zahl der gemessenen Datenpunkte
    'DECIM'   : numpy.int16, # DSP Dezimierungsfaktor
    'DSPFVS'  : numpy.int16, # DSP Version
    'GRPDLY'  : numpy.float64, #? DSP Gruppendelay
    'SI'      : numpy.int32, # Transformationsgrösse
    'BYTORDA' : numpy.int8, # Byteanordnung, 0 und 1
    'DTYPA'   : numpy.int8, # Zahlentyp, 0 = 32 Bit integer
    'AQ_mod'  : numpy.int8, # Aquirierungsmodus
    'DIGMOD'  : numpy.int8, # Filtertyp
    'NS'      : numpy.int64, # Anzahl der Scans
    'PULPROG' : str, # Name des Pulse-Programms
    'NUC1'    : str, # Beobachteter Kern
    'INSTRUM' : str, # Instrumentenname
    'WDW'     : numpy.int8, # Window function type
    'PH_mod'  : numpy.int8, # Phasing type: 0 = None, 1 = "Normal"
    'PHC0'    : numpy.float64, # Phase 0. Ordnung
    'PHC1'    : numpy.float64, # Phase 1. Ordnung
    'SSB'     : numpy.int64, # Sine bell shift
    'MC2'     : numpy.int64, # F1 Detektionsmodus (?)
    'EXP'     : str, # Name des Experiments
}

searchProc = {
    'SF'      : numpy.float64, # Im Prinzip sollte hier die Frequenz von 0ppm stehen.
    'FTSIZE'  : numpy.int32, # Beinhaltet offenbar ebenfalls die Zahl der Zahlen in fid - aber scheint zu stimmen.
    'SI'      : numpy.int32, # Transformationsgrösse, scheint auch in proc zu sein..
    'PHC0'    : numpy.float64, # Phase 0. Ordnung
    'PHC1'    : numpy.float64, # Phase 1. Ordnung
    'BYTORDP' : numpy.int8, # Byte-Ordnung der prozessierten Daten
    'DTYPP'   : numpy.int8, # Zahlentyp, 0 = 32 Bit integer
    'WDW'     : numpy.int8, # Window function type
    'SW_p'    : numpy.float64, # Spektrumbreite für proc
}


# Information about the bruker files are taken from the reference material 
#   and fileutil.cs from SpinWorks.
#
#
# Thanks to the authors of this great program and for providing those files
#   which were a great help to understand those files.

class BrukerNMR_1D(nmr_prototype.NMR_1D):
    sPath = ''
    
    acqu = {}
    proc = {}
    
    aq_mod = None # Aquirierungsmodus (s. acqu[AQ_mod]), wie er intern verwendet wird.
    
    AQ_mod_SEQ = 1
    AQ_mod_SIM = 2
    AQ_mod_DSP = 3
    
    f_custom_l_border = None
    f_custom_r_border = None
    
    y_values = {}
    
    o_workspace = None
    
    def __init__(self, path):
        super(BrukerNMR_1D, self).__init__(path)
        self._prepare()
        super(BrukerNMR_1D, self).load()
        
    def _prepare(self):
        self.o_workspace = lib.Workspace._the_instance
        
        self.readAcquisition()
        self.readProcessing()
        
        if self.getCustomLeftBorder() == None:
            self.setCustomLeftBorder(self.getLeftBorder())
            
        if self.getCustomRightBorder() == None:
            self.setCustomRightBorder(self.getRightBorder())
        
    def __getstate__(self):
        odict = {
            'sPath' : self.sPath,
            #'acqu' : self.acqu,
            #'proc' : self.proc,
            'aq_mod' : self.aq_mod
        }
        
        return odict
    
    def __setstate__(self, ndict):
        self.__dict__.update(ndict)
        
        # Re-Read.. While developmenent this is important!
        self._prepare()
        
    def report(self):
        print self.sPath
        print self.acqu
        print self.proc
        print self.aq_mod
        
    def readAcquisition(self):
        # Öffnet die acqus-Datei. Diese beinhaltet die für fid verwendeten Parameter.
        # acqu beinhaltet i.d.R. die gleichen, es sei denn, das Experiment wurde neu gestartet.
        filename = os.path.join(self.o_workspace.getNMRDir(), self.sPath, 'acqus')
        for line in open(filename, 'r'):
            if len(line) > 3:
                lineparts = list(line.partition('='))
                lineparts[0] = lineparts[0].strip().strip(' ><').lstrip('#$')
                lineparts[2] = lineparts[2].strip().strip(' ><')
                # Ende der Datei
                if lineparts[0] == 'END':
                    break    
                # Durchschauen der vorbereiteten Parameter
                elif lineparts[0] in searchAcqu:
                    self.acqu[lineparts[0]] = searchAcqu[lineparts[0]](lineparts[2])
        
        # Umdefinieren bzw Interpretieren der Bruker-AQ-Modi
        if self.acqu['AQ_mod'] == 1 or self.acqu['AQ_mod'] == 3:
            self.aq_mod = self.AQ_mod_SIM
        elif self.acqu['AQ_mod'] == 2:
            self.aq_mod = self.AQ_mod_SEQ
            
        if self.acqu['DECIM'] > 1:
            self.sq_mod = self.AQ_mod_DSP
            
    def readProcessing(self):
        # Öffnet die acqus-Datei. Diese beinhaltet die für fid verwendeten Parameter.
        # acqu beinhaltet i.d.R. die gleichen, es sei denn, das Experiment wurde neu gestartet.
        filename = os.path.join(self.o_workspace.getNMRDir(), self.sPath, 'pdata', '1', 'procs')
        for line in open(filename, 'r'):
            if len(line) > 3:
                lineparts = list(line.strip().partition('='))
                lineparts[0] = lineparts[0].strip(' <>').lstrip('#$')
                lineparts[2] = lineparts[2].strip(' <>')
                # Ende der Datei
                if lineparts[0] == 'END':
                    break
                # Durchschauen der vorbereiteten Parameter
                elif lineparts[0] in searchProc:
                    self.proc[lineparts[0]] = searchProc[lineparts[0]](lineparts[2])
                    
    def getExperiment(self):
        """ Return the experiment name """
        
        # HACK!
        # I don't know why -- but pythons seem to reset the name of the experiment
        # to "None" before this method get called. Just reread the values and we are
        # safe. Have to recheck this later.
        # /HACK!
        self._prepare()
        
        return self.acqu['EXP']
        
    def getOffset(self):
        """ Return the ppm-value of the first value. """
        return (self.acqu['SFO1'] / self.proc['SF'] - 1) * 1.0e6 + 0.5 * self.acqu['SW'] * (self.acqu['SFO1'] / self.proc['SF'])
        
    def getSpectralWidth(self, ppm = True):
        """ Return the Width of the Spectra.
        
        Keyword arguments:
        ppm -- True to return ppm, False to return Hertz
        
        """
        if ppm == True:
            return self.acqu['SW']
            
    def getLeftBorder(self):
        return self.getOffset()
        
    def getRightBorder(self):
        return self.getOffset() - self.getSpectralWidth()
        
    def getCustomLeftBorder(self):
        return self.f_custom_l_border
        
    def setCustomLeftBorder(self, border):
        self.f_custom_l_border = border
        
    def getCustomRightBorder(self):
        return self.f_custom_r_border
        
    def setCustomRightBorder(self, border):
        self.f_custom_r_border = border
        
    def getCustomSpectrumWidth(self):
        return abs(self.getCustomLeftBorder() - self.getCustomRightBorder())
        
    def getCustomSpectrum(self):
        return (
            self.y_values[self.ppm2i(self.getCustomLeftBorder()):self.ppm2i(self.getCustomRightBorder())], 
        )
        
    def i2ppm(self, i):
        """ Return the ppm value to a given point 
        
        Keyword arguments:
        i -- the Point to which the ppm value should be returned
        
        """
        resolution = len(self.y_values) / self.acqu['SW']
        ppm = self.acqu['SW'] - (i / resolution) - (self.acqu['SW'] - self.getOffset())
        return ppm
        
    def ppm2i(self, ppm):
        """ Return the point number i to a given ppm
        
        Keyword arguments:
        ppm -- the ppm value to which a point number should be returned
        
        """
        resolution = len(self.y_values) / self.acqu['SW']
        i = resolution * (self.getOffset() - ppm)
        return int(round(i, 0))
        
    def ppm2Hz(self, ppm):
        return  self.proc['SF'] * (ppm * + 1)
    
    def i2Hz(self, i):
        ppm = self.i2ppm(i)
        return self.ppm2Hz(ppm)
    
    def doFFT(self):
        """ Do the Fourier Transformationen. 
        ... okay, okay. It doesn't work yet. Just load the final spectra...
        
        """
        self.loadSpectra()
        #self.loadFid()
        
    def loadSpectra(self):
        """ Load the spectra """
        # Loads pdata/1/1r
        dataFile = open(os.path.join(self.o_workspace.getNMRDir(), self.sPath, 'pdata', '1', '1r'), 'rb')
        bytesRead = 0
        byteLength = 4
        dataPoints = self.proc['SI']
        
        values = numpy.zeros((self.proc['SI']), dtype = numpy.float)
        
        # Array vorbereiten
        if self.proc['DTYPP'] == 0:
            formatBE = '>l'
            formatLE = '<l'
        else:
            raise Exception("Unknown datatype: DTYPP %s" % (self.proc['DTYPP']))
        
        while bytesRead < dataPoints:
            byte = dataFile.read(byteLength)
    
            if self.proc['BYTORDP'] == 1:
                # Big Endian
                val = struct.unpack(formatBE, byte)[0]
            else:
                # Little Endian
                val = struct.unpack(formatLE, byte)[0]
            
            values[bytesRead] = val
            bytesRead += 1
        
        dataFile.close()
        self.y_values = self.normalize(values)
        
        # Save standard derivation - for peakpicking algorithm (among other usage)
        self.y_values_std = self.y_values.std()
        self.y_values_sino = min((self.y_values[0:2048].std(), self.y_values[dataPoints-2048:dataPoints].std()))
        
    def normalize(self, vals):
        maxima = vals.max()
        #length = len(vals)
        #valsRet = numpy.zeros(length, dtype = numpy.float)
        #
        #for i in range(0, length):
        #    valsRet[i] = numpy.float(vals[i])/maxima * 100
        
        valsRet = vals / maxima * 100
            
        return valsRet
            
    def loadFid(self):
        """ Load the FID (Free induction decay) """
        fid = open(os.path.join(self.o_workspace.getNMRDir(), self.sPath, 'fid'))
        # values = array.array('l')
        values = []
        
        byteread = 0
        bytelength = 4
        
        while byteread < self.acqu['TD']:
            byteread = byteread + bytelength
            byte = fid.read(bytelength)
            
            if self.acqu['BYTORDA'] == 1:
                # Big Endian
                val = struct.unpack('>l', byte)
            else:
                # Little Endian
                val = struct.unpack('<l', byte)
                
            print val[0]
        
        fid.close()

    def getValues(self):
        """ Return an array of y-values """
        return (
            self.y_values, 
        )
